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Abstract — Incoherence between sparsity basis and sensing basis 
is an essential concept for compressive sampling. In this context, 
we advocate a coherence-driven optimization procedure for vari- 
able density sampling. The associated minimization problem is 
solved by use of convex optimization algorithms. We also propose 
a refinement of our technique when prior information is available 
on the signal support in the sparsity basis. The effectiveness of the 
method is confirmed by numerical experiments. Our results also 
provide a theoretical underpinning to state-of-the-art variable 
density Fourier sampling procedures used in MRI. 

Index Terms — compressed sensing, variable density sampling, 
magnetic resonance imaging. 

I. Introduction 

Compressed sensing demonstrates that sparse signals can be 
sampled through linear and non-adaptive measurements at a 
sub-Nyquist rate, and still accurately recovered by means of 
non-linear iterative algorithms. The theory requires incoher- 
ence between the sensing and sparsity bases and a lot of work 
has thus been dedicated to design such sensing systems JT). 

In the present work, we concentrate on s-sparse digi- 
tal signals a. = (&i)-t <i<N G in an orthonormal 
basis U' = (tpi, i/)n) G <C NxN . The vector a con- 
tains s non-zero entries and its support is defined as S = 
{i : \cti\ > 0, 1 < i ^ N}. We denote a s G C s the vector 
made of the s non-zero entries of a. This signal is probed 
by projection onto m vectors of another orthonormal basis 
<t> = (<p u ...,<p N ) G C NxN . The indices of the selected 
vectors are denoted fi = {lx, . . . , l m } and 0^ is the m x N 
matrix made of the selected rows of <t>1\ where the symbol -t 
stands for the conjugate transpose operation. The measurement 
vector y G C m thus reads as 



= *3> G 

G C NxN . Finally, we aim at 
i-minimizatiorQ problem 

Ana. (2) 



y = An a with A n 

We also denote A = <t>tv|/ 
recovering a by solving the 

a = argmin ||ct||i subject to y 

aeC N 



imxN 



(1) 
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In this setting, common strategies focus on uniform random 
selection of the indices ij., . . . , l m . For signals sparse in the 
Dirac basis, a uniform random selection of Fourier basis 
vectors represents the best sampling strategy. Indeed, the Dirac 
and Fourier basis are optimally incoherent. Natural signals 
are however rather sparse in multi-scale bases, e.g. wavelet 
bases, not optimally incoherent with the Fourier basis. Many 
measurements are thus needed to reconstruct such signals 
accurately. This is for example the case in magnetic resonance 
imaging (MRI). To reduce the number of measurements, the 
authors in |3) rely on the fact that the energy of MRI signals is 
essentially concentrated at low frequencies. They thus propose 
to select Fourier basis vectors according to a variable density 
sampling profile selecting more low frequencies than high fre- 
quencies. This approach was shown to drastically enhance the 
quality of the reconstructed signals. This method is however 
essentially empirical and the reconstruction quality depends on 
the shape of the sampling profile used. Let us also mention 
that a line of justification for VDS was proposed in terms of 
the variable sparsity of the signals of interest as a function of 
scale in a wavelet sparsity basis JT], |4j. 

In this letter, we study VDS in the theoretical framework of 
compressed sensing. In Section fill we describe the latest com- 
pressed sensing results for sparse signals probed in bounded 
orthonormal system, and explain how they encompass variable 
density sampling procedures. In Section [Till we introduce a 
minimization problem for the coherence between the sparsity 
and sensing bases, whose solution provides an optimized 
sampling profile. This minimization problem is solved with 
the use of convex optimization algorithms. We also propose a 
further refinement of our technique when prior information is 
available on the signal support S. In Section ITVl we illustrate 
the effectiveness of the method through numerical simulations. 
We also provide a comparison of the Fourier VDS profile in 
the presence of prior information and corresponding recon- 
struction qualities, with the state-of-the-art VDS approaches 
used in MRI. Finally, we conclude in Section IVl 

II. Variable density sampling 

In the setting presented in Section [I] the compressed sensing 
theory demonstrates that if the sampling indices l\,...,l m are 
chosen randomly and independently according to a discrete 
probability measure P defined on {1, . . . , N}, then a small 
number m <C N of random measurements are sufficient for 
an exact reconstruction of a 0. 

Theorem 1 (Theorem 4.2, |0). Let A = <t>t\|/ e £Nxn^ md 

a G C N be a s-sparse vector such thaQ sgn (as) G C s is a 
random Steinhaus sequence. Assume that the sampling indices 



2 sgn (ag) £ C is the s-dimensional vector with entries onl \cti\, Vi £ 5". 
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= {7i,...,/ m } are selected randomly and independently 
according to a discrete probability measure P defined on 
{1, . . . , N}. Let y = A^a e C m and define 



i 



max 



1(^,^)1 



AT 1 / 2 l^i.KiV pV 2 (i) ' 
For a universal constant C > 0, if 

m > CjV/i 2 (P)s log 2 (6jV/e), 



(3) 



(4) 



f/ien a « the unique minimizer of the l\-minimization problem 
(0 with probability at least 1 — e. 

In the above theorem, the parameter p(P) stands for the 
mutual coherence between the measurement basis <t> and 
the sparsity basis U'. This value depends on the probability 
measure P and statisfies fJ,(P) iV -1 / 2 Q. The smaller 
the mutual coherence the smaller the required number of 
measurements for exact recovery. 

Let us highlight that with the selection procedure described 
in Theorem [T| the number of measurements is exactly m 
but one measurement vector might be selected more than 
once. This characteristic is not always suitable in practical 
applications, such as MRI, particularly in a VDS configuration. 
Indeed, a sensing basis vector <pi, whose associated probability 
of selection P(i) is high, will be selected multiple times thus 
reducing the quantity of information probed. To avoid this 
phenomenon, we propose another selection process. 

In the remainder, the sampling indices are selected accord- 
ing to an admissible sampling profile for m measurements. 

Definition 1 (Admissible sampling profile). A vector p = 



piV 



number m of measurements if pj £ (0, 1] for all 1 ^ j ^ N, 
and ||p||i = m. The set of all admissible sampling profiles for 
a number m of measurements is denoted V(m). 

Let p G V(m) be an admissible sampling profile, the 
sampling indices are selected by generating a sequence 
(5i, ...,Sn) 6 R w of independent Bernouilli random variables 
taking value or 1 and such that <5,; is equal to 1 with 
probability pi for all 1 ^ i ^ N. The set of selected indices is 
then defined as O = {I : Si = 1}. With the proposed sampling 
strategy, one measurement vector can be selected only once. 
The constraint that |p||i = m imposes that the number of 
measurements is m on average over realizations of a sequence 
(Si, (5 at). Note that for N 1, the variability of the number 
of measurements is negligible. 

As suggested in [2j, one can actually show that the recovery 
condition (O still holds with the coherence 



<*>-(*) 



1/2 



max 



\(<f>i,i>i)\ 



p 



1/2 



(5) 



The required elements of proof are provided in Appendix lAl 

III. Sampling profile optimization 

Let us assume that the number of measurements m is fixed. 
In order to recover the highest sparsity s possible, Theorem 
Q] shows that we should use the sampling profile p £ V(m) 
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Fig. 1, Top panels: probability of recovery e of s-sparse signals in the Haar 
wavelet basis (JV = 1024) as a function of the number of measurements m in 
the Fourier basis (left panel) and in the Hadamard basis (right panel). The dark 
dashed, dark continuous, and light continuous curves show the probability of 
recovery with a uniform sampling, an optimized variable density sampling, 
and the spread spectrum technique respectively. Curves on the left correspond 
to s = 50 and those on the right to s = 200. Bottom panels: light curves 
show the optimized sampling profile for m = 300 obtained with a sampling 
in the Fourier basis (left panel) and in the Hadamard basis (right panel). Dark 
curves show the values max^j^jv \ {<t>ii V'j)! 2 Ior all 1 ^ i ^ JV. 

minimizing the mutual coherence p(p). Therefore, we propose 
to solve the following optimization problem 



is an admissible sampling profile for a = 



arg mm 

(p,q)eK Nx2 



|B<z||oo + A||p-g-l|| 2 s.t. pelC T , (6) 



where A e [0,+oc), r e (0,1], IC T = {p e [t,1} N : ||p||i < 
m}, 1 £ M. N is the vector with all its entries equals to 1, 
p ■ q is the entry-by-entry multiplication between the vector 
p and q, \\ ■ H2 and || • are respectively the ^-norm and 
^oo-norrrH, and B £ C NxN is the diagonal matrix with entries 
maxi^j^jv I (4>ii Tpj) | 2 on tne diagonal, 1 ^ i ^ N. 

In the above problem, the term ||p ■ q — 1||| ensures that 
Pi ~ 1/qi for all 1 ^ i ^ N. The higher the value of the 
parameter A the further this constraint is enforced. In the limit 
where p ■ q = 1, we have HBqUoo = p 2 (p) confirming that 
problem (|6]l seeks to minimize the mutual coherence. Note 
that the minimization problem imposes that p belongs to the 
set K, T which is different from the set V(m). Consequently, 
we do not have necessarily ||p||i = m. However, we note that 
in practice the constraint ||p||i ^ m is always saturated for 
high enough values of A. 

To solve problem (O, we adopt the following procedure: 



Set t = and p(°) 
repeat 

q W argmin^Riv ||Bg 
p(* +1 ) <- argmin peRJ v ||p • q 
t<-t + l; 
until convergence 



+ X\\pW-q-l\\l; 
^ - 111 2 , s.t. p g K T ; 



Ei 



\xj I and \\x\\ 



■■ maxi^j^jv Fj I 
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Subproblems at step [3j and |4]i are convex problems. The 
subproblem at step [3]) is solved iteratively using a forward- 
backward algorithm and the one at stepHJi thanks to a parallel 
proximal algorithm Q. Both algorithms require the computa- 
tion of simple proximity operators. The computation of the one 
corresponding to ||B ■ essentially reduces to a projection 
onto an £i-ball (see Appendix IB1. This projection, as well as 
the one onto the ^i-ball of radius m, can be computed using 
the methocQ presented in j6). Note that for both subproblems, 
the computational complexity at each iteration is essentially 
driven by these projections for which the method in [6| has a 
worst-case complexity of 0(N log N). For N = 1024, as in 
the forthcoming experiments, the overall algorithm converges 
in at most a few seconds. Our procedure therefore easily scales 
to larger N. 

When prior information is available on the signal support S, 
we can refine our technique to find a sampling profile adapted 
to this support. Indeed, if the signal support S is known in 
advance then Theorem Q] applies with the coherence 

/ m ( m gjesK^V^ V 2 ,~ 
u(p, S) = — — max — . (7) 

We let the reader refer to Appendix [A] and equation © for 
more details. An optimized sampling profile associated with 
the set S can thus be obtained by substituting the diagonal 
matrix C 6 M. NxN with entries s~ 1 J2jes Iv&iV'j)! on me 
diagonal, 1 ^ i ^ N, for the matrix B in problem ©. 

IV. Experiments 

In order to evaluate the proposed method in a general 
setting, we conduct two experiments. For the first one, we 
choose the Haar wavelet basis as the sparsity basis U' and the 
Fourier basis as the sensing basis 0. We generate complex 
s-sparse signals of size N = 1024 with s e {50,200}. The 
positions of the non-zero coefficients are chosen uniformly at 
random in {1, . . . , N}, their phases are set by generating a 
Steinhaus sequence, and their amplitudes follows a uniform 
distribution over [0, 1]. The signals are then probed according 
to relation (fl~|i and reconstructed from different number of 
measurements m by solving the ^-minimization problem (f2]i 
with the SPGL1 toolbox [6]. For each value of to, the selected 
sensing basis vectors are chosen using the method described in 
Section [TT] using either a uniform density profile or the profile 
p obtained by solving problem (O with A = 0.05. Each time, 
the probability of recoverjQ is computed over 200 simulations. 
For the second experiment, the same setting is used but with 
the Hadamard basis as the sensing basis 4>. 

In order to evaluate our method when prior information is 
available on the support S, we perform a simplified MRI 
experiment. In this perspective, an in vivo brain image of 
size 256 x 256 was acquired on a 7 Tesla scanner (Siemens, 
Erlangen, Germany). As suggested in J5), we consider a 
Daubechies-4 wavelet basis as sparsity basis and decompose 

4 Code available at http://www.cs.ubc.ca/labs/scl/spgll 
5 Perfect recovery is considered if the l2-norm between the original signal 
x and the reconstructed signal x* satisfies: \\x — x*\\'2 ^ 10 J 1 1 1 1 2 . 




Fig. 2. Left panel: probability of recovery e of a MRI signal as a function 
of the number of measurements m. The light dot-dashed, light continuous, 
dark dot-dashed, dark continuous and dark dashed curves show, respectively, 
the probability of recovery obtained with: a uniform profile (a); the optimized 
profile obtained with the matrix B (b); the sampling profile proposed in (4) 
(e); the optimized profile obtained with the matrix C (c); the typical MRI 
profile used in (3) (d). Right panel: Optimized sampling profile obtained with 
the matrix C for m = 70 (dark continuous curve) in comparison with the 
one used in MRI (light continuous curve) and the one proposed in [4j (light 
dashed curve). 

each linJH of the brain image into this basis. The resulting 
vectors are then hard-thresholded at s = 50. All vectors but 
one are seen as a data set providing prior information on the 
support S of typical MRI signals. The average of the values 

YljeS I (^i' ' f° r eacn * m {!'•••! N}, serves to create 
the matrix C. The remaining signal, not considered in the 
data set, is considered as the signal under scrutiny, probed 
according to relation (Q]i and reconstructed from different 
number of measurements m by solving the l\ -minimization 
problem. For each value of to, the selected Fourier basis 
vectors are chosen using the method described in Section [TT] 
with: a uniform density profile (a); the optimized sampling 
profile p obtained with A = 0.05 and the matrix B (b) or with 
the matrix C (c); a typical sampling profile used in MRI (d) 
J51; the sampling profile proposed in j4| (e). 

Figure Q] shows the probability of recovery e of s-sparse 
signals as a function of the number of measurements for the 
first two experiments. The probability of recovery obtained 
with the spread spectrum technique is also presented Q, 10. 
Note that this technique, as for random Gaussian matrices, 
was proved to be universal, i.e., the number of measurements 
for the recovery of sparse signals is reduced to its minimum 
independently of the sparsity basis. One can note that the 
probability of recovery with the optimized sampling is always 
better than with the uniform sampling. With a sampling in the 
Fourier basis, one can also note that the recovery becomes 
almost optimal. Indeed, the number of measurements needed 
to reach a probability 1 of recovery is almost the same with 
the spread spectrum technique and with an optimized profile. 
These results confirm our theoretical predictions and illustrate 
the efficiency of variable density sampling. 

As illustration, Figure [T] also shows optimized sampling 
profiles obtained for the two sensing bases and to = 300 as 
well as the corresponding values of the diagonal entries of 
the matrix B. One can note that the shapes of the sampling 

6 Lines without any signal (background) are withdrawn. After this operation, 
151 lines are left. 

7 Note that the intrinsic parameters the sampling profiles (d) and (e) are 
manually chosen to obtained the best reconstructions. 
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profiles are highly correlated to the values in the matrix B. 

Figure [2] shows the probability of recovery e of the MRI 
signal as a function of the number of measurements. One 
can note that with the uniform sampling (a), the signal is 
recovered with probability 1 only when m = N. The results 
are slightly improved with the optimized profile (b) obtained 
with the matrix B. The sampling profiles (c), (d), and (e) 
drastically enhance the performance. Our optimized profile (c) 
obtained with the matrix C performs better than the profile (e) 
and similarly to the profile (d) typically used in MRI. These 
results provide a theoretical underpinning to VDS procedures 
used in MRI. It also shows that the refinement proposed for 
our technique can drastically enhances the performance of 
compressed sensing in practical applications. 

For illustration, Figure [2] also shows the sampling profiles 
(c), (d), and (e) for m = 70. One can notice that the profiles 
(c) and (d) are very similar to each other. This explains again 
the effectiveness of VDS profiles commonly used in MRI. 

V. Conclusion 

In the aim of optimizing variable density sampling profiles 
in the context of compressed sensing, we have introduced a 
minimization problem for the coherence between the sparsity 
and sensing bases. This problem is solved with the use of con- 
vex optimization algorithms. We have also discussed a refine- 
ment of our technique when prior information is available on 
the signal support in the sparsity basis. The effectiveness of the 
method is confirmed by numerical experiments. In particular, 
for signals sparse in a wavelet basis and probed in the Fourier 
domain, simulations show that our technique leads to optimal 
recovery. Indeed our technique gives similar probabilities of 
recovery as the spread spectrum method recently proved to be 
optimal. Our results also provide a theoretical underpinning to 
VDS procedures used in MRI. 

Appendix A 

The proof of this theorem follows exactly the method used 
to prove Theorem 4.2. in J2). The only difference resides in 
the estimate of the singular values of the operator Aq^Aqs 
(see Theorem 7.3, 0), where Afis E C mxs is the restriction 
of the matrix Aq to the columns indexed by S. 

~*NxN 



P = {Pi}l^N e V ( m )> 
NxN the diagonal matrix with 



Lemma 1. Let A = <t>t\|/ £ 

5 £ (0, 1/2], and define P £ 

1 /2 

entries p { on the diagonal, 1 ^ i ^ N. Assume that the m 
measurement vectors are selected according to p and suppose 
that s 2. For a universal constant C > 0, the normalized 
matrix A = P^ 1 A satisfies HAj^Ans — 1 1 
at least 1 - 2 3 / 4 s cxp f m&2 



^ 5 with probability 



CN u. 2 (p)s 



Proof: Let us denote Y = Aq S Aqs 
I G C sxs where a; G C lxs is the ■ 
proof starts by noticing that E [Y] = 
T,ti a W> - 1 = A l A s - I = 0. We 



EN r - T - 
i=1 6ia\ai- 

row of Ag. The 

n=lPi a i a ^ - 1 = 

can thus continue 



with the use of a symmetrization technique to bound the 
expected value of the norm of Y (see Lemma 6.7 in 0, 
or proof of Theorem 3.1 in JT]). Let (ei, ...,ejv) be a 



Rademacher sequence independent of (5i, ...,5n) and p ^ 2, 
then E||Y||" < 2*E|| Ei^n £iha\ai\\ p . Noticing that A ns 
has at most rank s, using Fubini's theorem, Rudelson's lemma 
(see Lemma 6.18, 0) conditional on (6\, . . . ,6n), and the 
Cauchy Schwarz inequality yields 



EIIYP < 2 3 / 4+p s (?■ 



p/2 



lA^Ans!!"! 



max ||a 4 ||2 P 

l\Oi— 1 



The previous equation is identical to equation (7.6) in the 
proof of Theorem 7.3 in 0. We can follow the same remain- 
ing steps of this proof to terminate ours. We still need however 
to provide a bound on maxi : 5 i= i ||a,|||. If the support S is 
fixed and known in advance, we have 



max <ij 9 max 

i:<5i=l l<i<N 



Pi 



m 



In the general case where S is unknown, we can write 



max a 



'112 



s max 



sN 2 

= — Kp) ■ 

m 



(8) 



(9) 



Appendix B 
The proximity operator of 7 1 1 B • || c 



7 > 0, is the unique 



solution of prQJCyMB.il (?) = argi™ij;gi» l/2||qr — x\\% + 

7ll Ba; l|oo- 

Proposition 1. For any q G R N , with B G R NxN defined 
as above, we have proJCyiig.ii (q) = q — jproj c (q/j), 
where C = {x eR N : || B _1 a;||i < 1} and proj c denotes the 
projection onto the set C. 

Proof: From Theorem 14.3 in |5|. we have q = 
P rox 7l|B.|U (?) +7P r ox 7 _i|| B .||^ (q/7), for all q G R N . In 
the previous relation, ||B • denotes the Fenchel conjugate of 
|| B • I Ioq. As B is a bijection (it is a diagonal matrix with strictly 
positive entries), one can show that (||B • ||oo)* = L c(') where 
be denotes the indicator function of the set C (Proposition 
13.20, 0). Finally, we have prox^-x y B . y ^ = prox 7 _i tc( . ) = 
proj c . Combining the last result with the first relation termi- 
nates the proof. ■ 
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